##### Step 1: load data and packages #####



library(tidyverse)
library(lubridate)
library(data.table)
library(haven)
library(fixest)
library(magrittr)
library(ggthemes)

base_dir <- "..." # Manually change this home directory
data_dir <- paste0(base_dir, "data/")

setwd(base_dir)

options(scipen=999)

#### table setup ####
setFixest_etable(digits = 3,
                fitstat = ~ n + r2)

setFixest_dict(c(rating="Fox News Rating, 4/15/09-10/15/09",
                 diff_rating="Change in Fox News Rating",
                 pre_rating = "Previous 6 mo. Rating",
                 attend09_max_pct = "\\% Tea Party Rally Attendance",
                 attend09_pct = "\\% Tea Party Rally Attendance",
                 attend09 = "Tea Party Rally Attendance (avg)",
                 attend09_max = "Tea Party Rally Attendance (max)",
                 praindum="Rain Prob. Decile",
                 popdum="Population Decile",
                 regiondum="Region",
                 rainrally = "Rain on 4/15",
                 pos_lin_fnc="FNC Channel Pos.",
                 pos_lin_msnbc = "MSNBC Channel Pos.",
                 fips = "County",
                 dmacode = "Designated Market Area (DMA)",
                 stid = "Cable system",
                 pop = "County Population"
                )
    )

tablestyle  = style.tex(main="aer",
                 fixef.suffix = " FEs",
                 fixef.where ="var",
                 fixef.title = "\\midrule",
                 stats.title = "\\midrule",
                 tablefoot=F,
                 yesNo="\\checkmark")


## MADESTAM ET AL CONTROLS
demo <- c("lmincome","unempl","lpopdens","rural","whitep","blackp","hispp","foreign","diff_unempl0509")
pastvotes08 <- c("rep2008","repvotes08","demvoted08","totvoted08","obama")
pastvotes06 <- c("rep2006","repvotes06","demvoted06","totvoted06")
pastvotes08pct <- c("rep2008","repvotes08_pct","demvoted08_pct","totvoted08_pct","obama")
pastvotes06pct <- c("rep2006","repvotes06_pct","demvoted06_pct","totvoted06_pct")
prains <- paste0("prain", 2:10) 
popdums <- paste0("popdum", 2:10) 
regionds <- paste0("regiond", 2:4)


##### APRIL 2009 PROTEST DATA FROM MADESTAM ET AL #####
protest <- read_dta(paste0(data_dir, "TPcounty_FS&RF.dta")) %>%
  as.data.table
protest[, fips := as.integer(fips)]

### add channel positions and controls ###

chpos <- readRDS(paste0(data_dir, "cable_county.rds"))

# cable system controls
cable_controls <- c("has_Fox.News.Channel", "has_MSNBC", "has_both", "has_neither", "has_fnconly", "has_msnonly", "num_channels", "num_bc")

chpos <- chpos[position_year == 2009, ]

chpos[,`:=` (has_both = as.numeric(chan_config==4),
             has_fnconly = as.numeric(chan_config==3),
             has_msnonly = as.numeric(chan_config==2),
             has_neither = as.numeric(chan_config==1))]

protest_cable <- chpos[protest, on = .(fips)]
protest_cable[,pos_diff := pos_lin_fnc - pos_lin_msnbc]

protest_cable[is.na(pos_lin_fnc), pos_lin_fnc := 0]
protest_cable[is.na(pos_lin_msnbc), pos_lin_msnbc := 0]

fox_ratings <- readRDS(paste0(data_dir, "nielsen_ratings.rds"))

ratings_2010 <- readRDS(paste0(data_dir, "nielsen_ratings_2010_primary.rds")) %>%
  .[channel == "FNC" & intab > 0]





##### Step 2: regressions of protests on channel positions #####



protest_models <- feols(attend09_pct ~ pos_lin_fnc + pos_lin_msnbc + rainrally + has_both + has_fnconly + has_msnonly +
  prain.[2:10] + popdum.[2:10] + regiond.[2:4] +
  csw0(lmincome + unempl + lpopdens + rural + whitep + blackp + hispp + foreign + diff_unempl0509,
       rep2008 + repvotes08_pct + demvoted08_pct + totvoted08_pct + obama + rep2006 + repvotes06_pct + demvoted06_pct + totvoted06_pct),
  data = protest_cable, weights = ~ pop, cluster = ~statefips)

mean_y1 <- mean(protest_cable$attend09_pct[setdiff(1:nrow(protest_cable), -protest_models[[1]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y2 <- mean(protest_cable$attend09_pct[setdiff(1:nrow(protest_cable), -protest_models[[2]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(protest_cable$attend09_pct[setdiff(1:nrow(protest_cable), -protest_models[[3]]$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table 3 -----

etable(protest_models,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3))),
    # file = paste0(table_dir, "protests_on_fnc_position.tex"),
    title = "Effect of Exposure to Fox News On Tea Party Rally Size",
    digits=3,
    digits.stats=2,
    style.tex=tablestyle,
    drop = c("(Constant)",cable_controls, demo, pastvotes08pct, pastvotes06pct),
    label = "table:protest-pos",
    group=list("\\midrule\nCable System Controls"=cable_controls,
              "County Demographics"=demo,
               "2008 and 2006 Voting"=c(pastvotes08pct, pastvotes06pct),
               "\\midrule\nRain Prob. Decile FEs" = "prain",
               "Population Decile FEs" = "popdum",
               "Region FEs" = "region"),
    replace = TRUE
  ) %>% print()





##### Step 3: regressions of ratings on protests #####



## DMA info for time zone
fipsdma <- fread(paste0(data_dir, "DMA-zip.csv")) %>%
  setnames(c("fips", "countyname", "state", "dmacode", "dmaname", "zipcode")) %>%
  .[,.(fips, dmacode)] %>%
  unique %>%
  .[,timezone := substr(as.character(dmacode), 1,1)] %>%
  .[,timezone := case_when(timezone == "5" ~ "ETZ",
                           timezone == "6" ~ "CTZ",
                           timezone == "7" ~ "MTZ",
                           timezone == "8" ~ "PTZ"
                          )] %>%
  .[is.na(timezone), timezone:= "PTZ"]   # alaska

## descriptive plot of ratings by timezone, day part, week day
## note M-Thu has highest ratings, but pacific TZ has high ratings at 5pm
## whereas eastern TZ has highest ratings at 8PM
fox_ratings[,wday := lubridate::wday(date, label=T)]
fox_ratings[,wday := fct_relevel(wday, "Sun", after=Inf)]
fox_ratings[,daypart := recode_factor(daypart, `M-F 07:00 AM - 09:00 AM` = "7-9 AM",
                                        `M-F 09:00 AM - 04:00 PM` = "9AM - 4PM",
                                        `M-F 04:00 PM - 05:00 PM` = "4-5 PM",
                                        `M-F 05:00 PM - 06:00 PM` = "5-6 PM",
                                        `M-F 06:00 PM - 07:00 PM` = "6-7 PM",
                                        `M-Sa 07:00PM - 10:00 PM; Su 06:00 PM - 10:00 PM` = "7-10PM",
                                        `M-Su 08:00 PM - 11:00 PM` = "8-11 PM")]
fox_ratings[,post := as.numeric(date >= mdy("04/15/2009"))]
fox_ratings[,taxday := as.numeric(date == mdy("04/15/2009"))]
fox_ratings[,days_since := time_length(date - mdy("04/15/2009"), unit = "day")]
fox_ratings[,weeks_since := floor(time_length(date - mdy("04/15/2009"), unit = "week"))]

# determine whether the county is in auto-collected markets (obs every day)
# or in diary markets (obs only in sweeps periods)
fox_ratings[, n_by_fips := .N, by = .(fips)]
fox_ratings[, diary_dma := as.numeric(n_by_fips <= 624)]

fox_ratings <- fox_ratings[fipsdma, on = .(fips), nomatch=0]
fox_ratings[,timezone:=factor(timezone, levels=c("ETZ", "CTZ", "MTZ", "PTZ"))]

## drop obs with no one in the panel
fox_ratings <- fox_ratings[intab > 0]

# use pre-tax-day ratings to construct weights
rat_by_time <- fox_ratings[post == 0,
  .(avg_rating = weighted.mean(rating, intab, na.rm=T)), by = .(timezone, wday, daypart, diary_dma)]

# create composite ratings by FIPS / week, weighting by tz-specific distribution computed above
rat_by_time[, weight := avg_rating / sum(avg_rating), by = .(timezone, diary_dma)]

fox_ratings <- fox_ratings[rat_by_time, on = .(timezone, wday, daypart, diary_dma)]

fox_ratings_prepost <- fox_ratings[,.(rating = 100*weighted.mean(rating, weight, na.rm=T),
                                     share = 100*weighted.mean(share, weight, na.rm=T),
                                     intab = weighted.mean(intab, weight, na.rm=T)),
                                  by = .(dmacode, fips, post, timezone, diary_dma)]

rat_pro <- fox_ratings_prepost[protest_cable, on = .(fips), nomatch=0]

rat_pro[, pre_rating := rating[which(post == 0)], by = .(fips)]
rat_pro[, diff_rating := rating - pre_rating, by = .(fips)]



# Print Figure D.2.1 -----

print(ggplot(
  aes(
      x = rating,
      y = attend09_pct,
      size = intab,
      weight = intab),
      data = rat_pro[post == 0]
    ) +
  geom_point(alpha = 0.4) +
  xlim(0, 7.5) +
  theme_bw() +
  theme(
    axis.text = element_text(size=12),
    axis.title = element_text(size=15),
    legend.title = element_text(size=15),
    legend.text = element_text(size=12)
  )+
  labs(x = "FNC Primetime Rating, 10/15/08 - 4/14/09", y = "Tea Party Rally Attendance (percentage of population)", size="Nielsen HH") +
  geom_smooth(method="lm", show.legend=F))



### 2010 summer ratings (04/01/2010-07/31/2010)

ratings_2010[, rating := rating * 100]  #rescale

rat_pro_2010 <- ratings_2010[protest_cable, on = .(fips)] %>%
  .[fox_ratings_prepost[post == 0, .(fips, pre_rating = rating, diary_dma)], on = .(fips)]

rat_pro_models_2010 <- feols(rating ~ attend09_pct +
  prain.[2:10] + popdum.[2:10] + regiond.[2:4] +
    pos_lin_fnc + pos_lin_msnbc + has_both + has_fnconly + has_msnonly + diary_dma +
  csw0(pre_rating,
       lmincome + unempl + lpopdens + rural + whitep + blackp + hispp + foreign + diff_unempl0509, 
       rep2008 + repvotes08_pct + demvoted08_pct + totvoted08_pct + obama + rep2006 + repvotes06_pct + demvoted06_pct + totvoted06_pct),
  data = rat_pro_2010, weights = ~ pop, cluster = ~statefips)

rat_pro_models_iv_2010 <- feols(rating ~ 
  prain.[2:10] + popdum.[2:10] + regiond.[2:4] +
    pos_lin_fnc + pos_lin_msnbc + has_both + has_fnconly + has_msnonly + diary_dma +
  csw0(pre_rating,
       lmincome + unempl + lpopdens + rural + whitep + blackp + hispp + foreign + diff_unempl0509, 
       rep2008 + repvotes08_pct + demvoted08_pct + totvoted08_pct + obama + rep2006 + repvotes06_pct + demvoted06_pct + totvoted06_pct) |
  attend09_pct ~ rainrally,
  data = rat_pro_2010, weights = ~ pop, cluster = ~statefips)

setFixest_dict(c(rating="Fox News Rating, Summer 2010",
                 diff_rating="Change in Fox News Rating",
                 pre_rating = "Pre-April 2009 Rating",
                 attend09_max_pct = "\\% Tea Party Rally Attendance",
                 attend09_pct = "\\% Tea Party Rally Attendance",
                 attend09 = "Tea Party Rally Attendance (avg)",
                 attend09_max = "Tea Party Rally Attendance (max)",
                 praindum="Rain Prob. Decile",
                 popdum="Population Decile",
                 regiondum="Region",
                 rainrally = "Rain on 4/15",
                 pos_lin_fnc="FNC Channel Pos.",
                 pos_lin_msnbc = "MSNBC Channel Pos.",
                 fips = "County",
                 dmacode = "Designated Market Area (DMA)",
                 stid = "Cable system",
                 pop = "County Population",
                 diary_dma = "Nielsen Diary Market"
                )
    )



mean_y1 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_2010[[1]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y2 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_2010[[2]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_2010[[3]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y4 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_2010[[4]]$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table 4 -----

print(etable(rat_pro_models_2010,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4))),
    # file = paste0(table_dir, "ratings_2010_on_protests_ols.tex"),
    title = "Effects of Tea Party Rally Size on Subsequent FNC Ratings (OLS)",
    tex=TRUE,
    digits=3,
    digits.stats=2,
    style.tex=tablestyle,
    drop = "(Constant)",
    label = "table:ratings2010-ols",
    group=list("\\midrule\nCable System Controls"=c(cable_controls, "Nielsen Diary Market"),
              "Cable Positions" = "Channel Pos",
              "County Demographics"=demo,
               "2008 and 2006 Voting"=c(pastvotes08pct, pastvotes06pct),
               "\\midrule\nRain Prob. Decile FEs" = "prain",
               "Population Decile FEs" = "popdum",
               "Region FEs" = "region"),
    replace = TRUE
  ))



mean_y1 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_iv_2010[[1]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y2 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_iv_2010[[2]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_iv_2010[[3]]$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y4 <- mean(rat_pro_2010$rating[setdiff(1:nrow(rat_pro_2010), -rat_pro_models_iv_2010[[4]]$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table D.3.1 -----

etable(rat_pro_models_iv_2010,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4))),
    # file = paste0(table_dir, "ratings_2010_on_protests_iv.tex"),
    title = "Effects of Tea Party Rally Size on Subsequent FNC Ratings (2SLS)",
    digits=3,
    digits.stats=2,
    style.tex=tablestyle,
    drop = "(Constant)",
    label = "table:ratings2010-2sls",
    group=list("\\midrule\nCable System Controls"=c(cable_controls, "Nielsen Diary Market"),
              "Cable Positions" = "Channel Pos",
              "County Demographics"=demo,
               "2008 and 2006 Voting"=c(pastvotes08pct, pastvotes06pct),
               "\\midrule\nRain Prob. Decile FEs" = "prain",
               "Population Decile FEs" = "popdum",
               "Region FEs" = "region"),
    replace = T,
    placement = "H"
  ) %>% print()





##### Step 4: placebo test #####



rainfall_placebo <- feols(rainrally ~ pos_lin_fnc + pos_lin_msnbc + has_both + has_fnconly + has_msnonly +
                            prain.[2:10] + popdum.[2:10] + regiond.[2:4] +
                            csw0(lmincome + unempl + lpopdens + rural + whitep + blackp + hispp + foreign + diff_unempl0509,
                                 rep2008 + repvotes08_pct + demvoted08_pct + totvoted08_pct + obama + rep2006 + repvotes06_pct + demvoted06_pct + totvoted06_pct),
                          data = protest_cable, weights = ~ pop, cluster = ~statefips)

rainfall_placebo

# Print Table F.1 -----

etable(rainfall_placebo,
       # file = paste0(table_dir, "rainfall_on_fnc_position.tex"),
       title = "Placebo Test of Channel Positions On Rainfall on Tax Day, 2009",
       digits=3,
       digits.stats=2,
       style.tex=tablestyle,
       drop = c("(Constant)",cable_controls, demo, pastvotes08pct, pastvotes06pct),
       label = "table:rainfall-placebo",
       group=list("\\midrule\nCable System Controls"=cable_controls,
                  "County Demographics"=demo,
                  "2008 and 2006 Voting"=c(pastvotes08pct, pastvotes06pct),
                  "\\midrule\nRain Prob. Decile FEs" = "prain",
                  "Population Decile FEs" = "popdum",
                  "Region FEs" = "region"),
       replace = TRUE,
       placement = "H"
) %>% print()


